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& ABSTRACT 

A main science goal for the Large Synoptic Survey Telescope (LSST) is to measure the 
i— i cosmic shear signal from weak lensing to extreme accuracy. One difficulty, however. 

is that with the short exposure time (~15 seconds) proposed, the spatial variation of 
the Point Spread Function (PSF) shapes may be dominated by the atmosphere, in 
addition to optics errors. While optics errors mainly cause the PSF to vary on angular 
scales similar or larger than a single CCD sensor, the atmosphere generates stochastic 
structures on a wide range of angular scales. It thus becomes a challenge to infer the 
I multi-scale, complex atmospheric PSF patterns by interpolating the sparsely sampled 

stars in the field. In this paper we present a new method, PSFENT, for interpolating 
the PSF shape parameters, based on reconstructing underlying shape parameter maps 



with a multi-scale maximum entropy algorithm. We demonstrate, using images from 
the LSST Photon Simulator, the performance of our approach relative to a 5th-order 
• • polynomial fit (representing the current standard) and a simple boxcar smoothing 

technique. Quantitatively, PSFENT predicts more accurate PSF models in all scenarios 
and the residual PSF errors are spatially less correlated. This improvement in PSF 
interpolation leads to a factor of 3.5 lower systematic errors in the shear power spec- 
trum on scales smaller than ~ 13', compared to polynomial fitting. We estimate that 
with PSFENT and for stellar densities greater than ~l/arcmin 2 , the spurious shear 
correlation from PSF interpolation, after combining a complete 10- year dataset from 
LSST, is lower than the corresponding statistical uncertainties on the cosmic shear 
power spectrum, even under a conservative scenario. 

Key words: cosmology: observations - gravitational lensing - atmospheric effects - 
methods: data analysis - techniques: image processing - surveys: LSST 
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1 INTRODUCTION 

Gravitational lensing is the physical phenomenon where 
gravitational fields perturb the trajectory of light rays and 
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therefore distort observed images | |Schneider|p92| ). In par- 
ticular, the study of weak gravitational lensing involves mea- 
suring the statistical properties of an ensemble of distorted 
galaxy images (see e.g. Bartelmann & Schneider 20011. 



Weak gravitational lensing is, in principle, one of the most 
powerful probes of dark matter and dark energy. By measur- 
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ing, at different redshifts, the statistical distortion of back- 
ground galaxies due to large scale cosmic structures - the 
"cosmic shear" - it is possible to place extremely tight con- 
straints on the nature of dark energy ( [Jain &: Seljaklp97| 
Hu fc Tegmark|1999 1 



et al.|2010l, the PSF anisotropy in long exposure data can 



For observations to date, the accuracy of cosmic shear 
measurements has been mostly limited by the statistical 
variation of random galaxy shapes in the relatively small 
sky areas studied ( |Hoekstra et al. | |2006| |Semboloni et al" 
20061 IHetterscheidt et al.||2006| IJarvis et al.||2006| IHetter- 



scheidt et al. 2007; Benjamin et al. 2007; Schrabback et al 



20101. However, in future wide-field weak lensing surveys 
such as those plann ed with the Dark Energy SurveyJ^LSSl^] 



llvezic et al. 



20081, and Euclid]^] extremely large datasets 
will greatly reduce the statistical errors, making these ex- 
periments systematics-limited. 

A major source of systematic error in weak lensing 



comes from our incomplete knowledge of the PSF (Paulin- 
Henriksson et al.|2008 1. To account for the effect of the PSF 



on observed galaxy images, we need a model for it at every 
galaxy position; these models can be constrained by images 
of stars, which provide noisy estimates of the PSF shape 
that are more sparsely distributed than the galaxies. This is 



the "PSF interpolation problem". In an earlier paper ( Chang 
et al.|[2012 hereafter C12), we quantified the spatial varia- 



tion in the PSF shapes for a typical LSST 15-second expo- 
sure due to various physical effects, such as optics misalign- 
ments and atmospheric turbulence. We found that, although 
most of the PSF anisotropy due to instrumental effects varies 
smoothly over the field of view, the atmospheric turbulence 
can generate PSF spatial variation on a wide range of scales 
in these short exposures, with patterns that do not repeat 
over time. This poses a new PSF interpolation challenge 
quite different from that faced by previous studies, which re- 
lied on images with longer exposure times and/or contained 
large instrumental effects that dominate the errors. 

These atmospheric features have been observed in short 
exposure (~ 10 seconds) images by, for example, Wittman 
( |2005| ) and |Heymans et al.| ( |2012| hereafter H12); in partic- 
ular, H12 pointed out that such high frequency, turbulence- 
induced spatial PSF variations in single short exposures may 
lead to systematic errors in shear measurements at levels 
that cannot be ignored for future weak lensing surveys, were 
existing PSF interpolation techniques to be used. 

In weak lensing analyses to date, the scheme used most 
often to interpolate the PSF between sparsely sampled stars 
has been to fit a low order two-dimensional spatial polyno- 
mial function to the stars' shape parameters. |Van Waer-| 
beke et al. ( 2002 1 claimed that 2nd-order polynomials are 



sufficient to model the PSF anisotropy variation across a 
typical CCD sensor of size ~ 10', while other studies iden- 
tified possible drawbacks of simple polynomial fitting and 
more sophisticated models have been suggested (e.g.|Massey| 
et al.|2002|IHoekstra|20lMl|Van Waerbeke et aL|2005||Berge1 



et aLI 2012 1. In general, being dominated by instrumental 



effects, which primarily generate large-scale, smooth fea- 



tures (Jarvis & Jain 2004 Rhodes et al. 2007 Schrabback 
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be modelled reasonably well using low order interpolation 
methods, even though these patterns are only sparsely sam- 
pled by the relatively low density of stars in the field. How- 
ever, for future synoptic wide-field surveys, such as LSST, 
where high cadence imaging is required, short exposures are 
inevitable. It is therefore important to re-examine the tra- 
ditional PSF interpolation techniques, and develop new in- 
terpolation algorithms that are better suited for these data. 

A deliberate effort is being made to study in detail 
the PSF patterns for LSST by carrying out end-to-end, 
photon-by-photon image simulations using the LSST Pho- 
ton Simula tor (PHoSlM; |Peterson et al.|2012|[2009{|Connolly| 



et al 



2010| ).P|The ray-tracing procedure includes models for 



the instrument response, telescope optics, and, most impor- 
tantly for our purposes, the atmosphere. For a full descrip- 
tion of the atmospheric model and quantitative comparison 
of the model against real data, we refer the reader to |Peter^| 
|son et aL| ( |2012| hereafter P12). In this work, our focus is on 
studying the PSF interpolation problem, for which we re- 
quire simulated data containing realistic atmospheric PSFs, 
with a range of strengths and spatial scales. The PhoSim 
images meet these criteria, as we demonstrated in P12 by 
comparing the relevant PSF characteristics in our simulated 
images with those seen by H12 in short exposure images 
taken by the MegaCam wide-field camera (1 degree 2 ) on 
the Canada- France-Hawaii Telescope (CFHT). In this work 
we therefore test our new algorithm on PhoSim simulated 
images, for which we know the underlying PSF spatial vari- 
ation. 

Working with simulated images is vital for this par- 
ticular task, since existing data that are suitable for our 
tests are very limited - our tests require wide-field, short- 
exposure images with high stellar densities, and taken under 
a wide range of atmospheric conditions. Simulations, with 
sufficiently high fidelity, allow us to test the absolute, as 
well as relative, accuracy of our algorithm in a more con- 
trolled fashion and with higher statistics. Our PhoSim ap- 
proach can be seen as the next step beyond that taken by 



Berge et al. (20121, who re-sampled PSF patterns observed 
in Subaru images. Here, we use PhoSim to generate large 
numbers of predicted LSST images with the relevant expo- 
sure time, 15 seconds - effectively amplifying the data taken 
with CFHT for a more rigorous testing. 

The primary aim of this paper is to introduce a new 
PSF interpolation technique, and test it under controlled 
conditions. In both real and simulated PSF patterns we ob- 
serve structure due to the atmosphere on many different 
angular scales: this motivates us to model the underlying 
anisotropy maps using a range of different-sized smoothing 
kernels. We infer the pixel values of these underlying maps 
from the noisy, sparsely sampled stellar shape data given 
a non-committal entropic prior. We quantitatively compare 
this new maximum entropy method with two other meth- 
ods that represent the current standards, and investigate the 
performance of all three using high-fidelity simulations over 
a range of observing conditions. 

The structure of this paper is as follows. In Section [2] 
we review briefly the relevant weak lensing theory and the 
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PSF interpolation problem for weak lensing. We then de- 
scribe, in Section [3] our new PSF interpolation method and 
give arguments for why it is well-suited for the particular 
problem at hand. In Section [4] we use simulated images to 
quantify the performance of our new method against two 
strawman PSF interpolation techniques. We also define the 
metrics that we use to quantify the performance of a given 
interpolation method. The results of this programme are 
presented in Section[5] We then discuss their implications in 
terms of the systematic errors in cosmic shear measurement 
for future weak lensing surveys and make suggestions for 
further improvements on psfent in Section [6] We conclude 
in Section [7] 



2 WEAK LENSING AND PSF 
INTERPOLATION 

In the weak lensing regime, the effect of gravitational lens- 
ing is to add a small offset to the intrinsic ellipticity of 
each galaxy, where ellipticity is typically defined as a two- 
component complex "spinor", e 



ei + i£2 (see e.g. Schnei- 
der et al.|2002 \ . The resulting observed ellipticity is a noisy 
but unbiased estimator of the applied shear. We then con- 
structs certain statistics from these shear estimators to infer 
cosmology. One of the most popular statistics is the two- 
point ellipticity correlation function £±(0) for the ensemble 
of galaxies: 

U(0) = (et(0o)e t (0o + 6)) ± <£ x (0 o )e x (ft) + 9)) , (1) 

where the angle brackets indicate an average over galaxy 
pairs separated by (with one galaxy located at 8o) and the 
subscripts t and x indicate an isotropised decomposition of s 
along the line connecting a certain pair of galaxies. The shear 
correlation functions predicted from cosmology, compared 
with these observed galaxy ellipticity correlation functions, 
provide a route by which the cosmological parameters can 
be inferred. 

The major challenge in a ground-based weak lensing 
analysis is to account for the instrumental and atmospheric 
PSF contribution to the observed galaxy shapes, such that 
these effects do not systematically contaminate the shear 
signal one wishes to measure. This involves "deconvolving" 
(approximately) the PSF from the galaxy images, where 
the PSF at each galaxy's location is "interpolated" from the 
shapes of nearby stars. 

A wide range of algorithms have been developed to 
model the shapes of the galaxies and stars, and to perform 



PSF deconvolution in the noisy data (see e.g. 


Hey mans et al. 


2006[ |Massey et al.||2007[ |Bridle et al.||2010[ 


Kitching et al. 


2012a|b for a summary of the various methods). However, 



to date, the PSF interpolation problem has been taken to be 
of secondary importance, since the PSF ellipticity patterns 
in existing images appear to be largely instrumental in ori- 
gin, somewhat repeatable, and well-modelled by smoothly- 
varying functions such as low-order polynomials. However, 
as we showed in P12 and C12, this may not be the case 
for future instruments such as LSST, which are specifically 
designed for weak lensing and have extremely tight require- 
ments on the instrument-induced PSF anisotropy. In these 
circumstances, the atmospheric effects, which used to be 
subdominant to instrumental effects, now become one of the 



key components in determining the PSF shape and the PSF 
spatial variation. This implies that the PSF shapes may no 
longer be smoothly varying across the field and modelling 
the PSF variation with such assumptions may be problem- 
atic. In addition, since the atmospheric effects are more pro- 
nounced in short exposures, datasets such as LSST, which 
are composed of sets of multi-epoch short exposure images 
instead of one long exposure, may suffer more from the at- 
mospheric effects in single exposures when constructing the 
PSF model from interpolation. 

To illustrate this PSF interpolation challenge, we show 
in Figure [I] two single-component (ei), model-subtracted 
stellar ellipticity maps from 74-second exposure images 
taken with the CFHT MegaCarrj^] on a dense stellar field 
(~ 7 stars per arcmin 2 ). The two images used to construct 
Figure [I] were taken on the same patch of sky but in two 
different nights, which appear to have very different atmo- 
spheric conditions. The data in Figure [I] are taken from two 
of the catalogues used in H12, so we refer to their paper for 
further details of the dataset. Note that a 2nd-order polyno- 
mial was subtracted from the raw ellipticity measurements 
(as explained in H12) - this accounts for most of the instru- 
mental PSF contribution, but may also have removed some 
large scale atmospheric features. These images demonstrate 
that PSF spatial patters contain the characteristic high fre- 
quency structures from the atmosphere, which is the main 
motivation for our new PSF interpolation method. 



3 PSFENT: A MULTI-SCALE INFERENTIAL 
INTERPOLATION METHOD 

As seen in Figure [T] the atmospheric PSF anisotropy pat- 
terns can contain structure on a range of angular scales, with 
both patchy and striped features. This motivates us to look 
for flexible functions with which to model this spatial vari- 
ability, which we can then fit to the sparsely sampled stellar 
PSF shape data in any given situation. 



3.1 Interpolant model, and the likelihood function 

Throughout this paper, we test the performance of our new 
interpolation method by interpolating two shape parame- 
ters: ei and £2- In practice, a full weak lensing pipeline will 
require interpolation of several other shape parameters as 
well (e.g.PSF size, matrix elements of the shear polarisabil- 



ity in KSB-type methods (Kaiser et al. 19951 etc). We do 
not do a full analysis on these other parameters, but demon- 
strate in Appendix [B] that the our methodology is easily 
generalised to other parameters. 

For the atmosphere-induced PSF, the two components 
of the ellipticity, ei and £2, are expected to be independent 
from each other in each exposure, while both varying with 
similar amplitudes and spatial structures between different 
exposures. (We demonstrate this point with simulations in 
Appendix [c]) This is because the physical mechanism that 
generates these ellipticity values - the refraction by turbu- 
lent cells - does not have a preferred direction, but does have 
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Figure 1. Examples of the residual PSF ellipticity (ei) patterns in ~1 degree 2 sky regions, as observed in a dense stellar field with the 
CFHT MegaCam. The exposure time for these images is 74 seconds. As described in H12, a 2nd-order polynomial model has been fitted 
and subtracted from the raw stellar ellipticities to remove the optics contribution. 



a certain magnitude and spatial power spectrum dictated by 
the local atmospheric parameters ( de Vries et al.|2007 l. As 
a result, we treat the two components of complex PSF el- 
lipticity as independent fields with the same priors. These 
two fields, s\(x, y) and £2(x, y), are to be reconstructed from 
sparse, noisy, stellar shape data e° b | and e\ jf . Here, k runs 
from 1 to the number of stars observed, N s tar- Casting the 
PSF interpolation problem as an image restoration prob- 
lem in this way allows us to properly take into account the 
observational errors on the measured star shapes, and prop- 
agate those errors into uncertainties on the interpolant. An 
iterative likelihood fit is performed: at each step, the two el- 
lipticity components of each stellar image are predicted from 
the model underlying ellipticity fields and compared to the 
measured stellar ellipticities. 

Both the predicted and observed data are inputs to the 
likelihood function. Under the assumption of uncorrelated 
Gaussian stellar shape uncertainties crfc, this can be writ- 
ten for e.g.the first ellipticity component as the following 
probability distribution (PDF) for the data: 



(2k) 



" /3 n fc *i.* 



x exp 



obs 
e l,fc 



■ Ei(xk,yk;hi) 



(Tl.k 



(2) 



and likewise for £2. 

In Equation [2] hi is a parameter vector that represents 
the model. It is the components of this parameter vector 
(and its companion h,2 for E2) that we vary to fit the stellar 
shape data. We choose to parameterise the flexible interpo- 
lation functions e%[x, y) and £2(1, y) with pixelated grids on 
the sky. We compute the predicted ellipticity at the fcth star 
position, si(xk,yk',hi), by linear interpolation between the 
neighbouring pixels: we choose the pixel scale of each grid 
on our maps such that each pixel contains approximately 
1 target point, on average, such that the linear interpola- 
tion choice does not affect the final prediction. For our test 
data, we fix the model map sizes at 80 x 80 pixels for each 
ellipticity component. 



Such free-form discretised functions like h\ and I12 have 
as many parameters to be inferred as there are pixels in the 
grids. However, we would like to impose some smoothness 
on these maps, such that structure on a range of angular 
scales can be predicted. We do this by constructing each map 
from a weighted sum of seven "hidden" maps, each convolved 
with a Gaussian "Intrinsic Correlation Function" (ICF) of a 
different angular scale: the result is known as the "visible" 
map. This procedure provides an efficient way of introducing 
smooth, correlated structure on a variety of angular scales. 
In our notation, h stands for "hidden." We therefore have 
7 x 80 x 80 = 44, 800 hidden pixel values to vary during the 
fit, for each ellipticity component. The convolutions with the 
ICF kernels reduce the effective number of free parameters, 
but even so many of these will still turn out not to be con- 
strained by the few hundred data points in the field. The 
choice of prior PDF for the parameters in the h is therefore 
important. 



3.2 The entropic prior PDF 

We take the pixel values of the hidden images to be uncorre- 
lated by construction, and assign a positive-negative entropic 



prior for them (Maisinger et al. 2004 hereafter MHL04) 



This has the effect of suppressing structure in the maps 
unless it is required by the data. In this way we give the 
method plenty of flexibility to fit the data well, but reg- 
ularise to avoid over-fitting. Such a multi-scale maximum 



entropy method was first used by Weir (19921, and is im- 



plemented in the publicaly available MemSys4 code (Gull 
& Skilling 1999). We illustrate the construction of a multi- 
scale ellipticity component map in Figure [2] showing how a 
range of different features on different angular scales can be 
modelled. 

This model is similar in both essence and outcomes to 
one comprising a pixelated map and its "a trous" wavelet 
transform, as shown in some detail by MHL04. This wavelet 
transform can also be written as a set of convolutions; the 
implementation of MHL04 works well when making maps 
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Figure 2. Illustrating the psfent multi-scale ellipticity map 
model. "Hidden" maps (left column, greyscale) are convolved with 
Gaussian kernels (ICFs) of exponentially-decreasing size (middle 
column) to make seven component "visible" maps (right, orange). 
These component maps are then weighted and summed to make 
the final model ellipticity map, which can be interpolated linearly 
onto any target position to predict the PSF shape there. 



of the CMB temperature anisotropies, which also exhibit 
patchy features on a range of angular scales. Just like wavelet 
basis functions, the Gaussian ICFs we use have characteristic 
angular scales Wi that increase approximately exponentially: 
the smallest is a single map pixel (wo — 1), while the largest 
is approximately wg = 2 6 = 64 pixels in size. The ICFs are 
normalised to unit volume. 

The entropic prior PDF for a single hidden pixel value 
hi takes the following form 

PT(hi\mi) =exp[aS(hi)] , (3) 
S(hi) =j>i- 2m, - h t log f ^ 2 ~M , (4) 

where tpi — (hf + 4mf) 1/ ' 2 . This distribution peaks at zero 
and is symmetric. The "regularisation constant" a is a (nui- 
sance) hyperparameter that parametrises the prior distribu- 
tion, a is inferred from the data via the Bayesian evidence 
internally by MemSys4, and controls the final importance 
of the prior relative to the data. The "model" values rrii 
(which we take to be constant over each hidden image) are 
also hyper-parameters, that determine the ease with which 
structure develops at each resolution scale: the smaller the 
value of rrii, the stronger the suppression of features at that 
scale. 

3.3 Informing the prior 

At this point we might ask whether we can inject any more 
information into the problem by choosing values of the rrii 
to reflect the statistical properties of the simulated atmo- 
spheric PSF patterns. For a given pixel value, the entropic 
prior has approximate width m ~ <t^/2, where au is the rms 
width of an approximating Gaussian (MHL04). This sug- 
gests that a possible algorithm for assigning the prior width 
at a particular resolution scale is to consider the variance of 
pixel histograms of low noise "true" multi-scale ICF hidden 
ellipticity maps at that scale. 

To do this, we generate a special set of 100 simulated 
images with an ultra-high density of stars and low noise. 
The range of PSF patterns in these simulations is consis- 
tent with the simulations used in the main analyses and is 
described in more detail in Section [4.1| These unphysically 
dense star fields allow us to access the true PSF elliptic- 
ity pattern expected in the short exposures. The simula- 
tions are then run through psfent to be effectively "decom- 
posed" into the 7 different scales, corresponding to the 7 
ICFs. We generated one histogram for each of the 7 resolu- 
tion scales that contain the pixel values (ei and £2) in the 
hidden images for all 100 atmospheric PSF patterns on those 
scales. These histograms, as shown in Figure [3] are more 
peaked, and have broader wings, than a Gaussian distribu- 
tion. The entropic prior PDF (in blue), while still imperfect, 
is a slightly better approximation to these distributions. We 
found the rms widths of the average histograms at the dif- 
ferent scales increase from the smallest scale to the largest 
scale as: m, = [0.0003, 0.0002, 0.0027, 0.0050, 0.0064, 0.0089, 
0.0188]. There is very little power on the smallest two scales 
and nonlinear increase from the remaining middle to large 
scales. We adopt these numbers as the width for the en- 
tropic priors throughout the rest of the paper. We find that 
informing the priors in this way largely improves the accu- 
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racy of the PSF model constructed by psfent, compared 
to using flat, non-informative priors. We also find that us- 
ing the default prior set in MemSys4 (which was optimised 
for CMB temperature map reconstruction) results in well- 
behaved PSF models, although quantitatively worse than 
the priors constructed from realistic simulations described 
abov(0 

The standard deviation of the width of the grey his- 
tograms is approximately 30 - 40 % the mean width of the 
same histograms. This variation is due the variation in the 
atmospheric conditions in our simulations. For further opti- 
misation of psfent, we could also use a narrower range of 
weather conditions and carry out the same process multiple 
times to derive priors as a function of different atmospheric 
parameters such as seeing and wind speed. 

As an aside, we note that the process of assigning prior 
widths for the different angular scales plays a very similar 
role to the assignment of the range parameter in the co- 
variance function of the Gaussian process at the heart of a 
Kriging interpolation ( |Berge et al.|2012 1 . The Kriging range 
parameter could also be derived by inspecting large numbers 
of high density simulated starfields, to capture the informa- 
tion present in the data-constrained atmospheric turbulence 
model in its most useful form. 



3.4 Estimating the posterior PDF 

The posterior PDF, e.g.Pi(hj |{e° bs }), for the hidden pixel 
values given the data can be approximated by a multivari- 
ate Gaussian, centred at the maximum posterior point. The 
maximum posterior maps provide our best estimates for the 
PSF shape parameters at any target point in the field. Sam- 
ple maps can be drawn from the posterior PDF in order to 
provide approximate uncertainties on these estimated PSF 
parameters. We find that 100 sample maps provide a suf- 
ficiently accurate standard deviation map, which we use 
for the uncertainties on the predicted PSF ellipticity esti- 
mates. This Gaussian approximation (including the maxi- 
mum of the posterior, the covariance matrix of the param- 
eters, and the associated evidence) is computed using the 
MemSys4 code, available on request from MaxEnt Data 
C onsultants^] Det ails of the implementation can be found 
in|Gull fc Skilling|(p99|). 



4 SIMULATION AND ANALYSIS 
4.1 Simulations 

As discussed in Section [I] since the existing datasets are in- 
sufficient for us to perform a systematic test with psfent, we 
depend on PhoSim to generate a large number of simulated 



images. We refer to P12, Peterson et al. ( 2009 1 and Connolly 



et al. (20101 for a complete description of PhoSim and in 



e We calculate the <t[£psf] and 5"g ya pgF values in the case of 
using psfent with the MemSys4 priors to interpolate the PSF 
ellipticities for the same set of simulations in Section [5.1| (See Sec- 
tion |4.3| for the definition of these two statistics.) Quantitatively, 
using the MemSys4 priors increases (rfepgp] by 3% and o-^ ya pgF 
by 29% compared to using the priors derived from simulations. 
7 http://www.maxent.co.uk 



specific the details of the atmospheric model. Here, to facili- 
tate comparison with parallel PSF interpolation studies, we 
provide a very brief overview of the PhoSim atmospheric 
model. 

At the heart of the PhoSim atmospheric model is a sys- 



tem of seven-layer frozen Kolmogorov screens ( Kolmogorov 
1992"l|Lane, R. G.|1992| ). These screens are constructed with 



density fluctuations obeying a full three-dimensional Kol- 
mogorov energy spectrum of E(k) oc fc~ 5//3 , with values 
for the mean seeing, inner and outer turbulence scale as- 
signed to each individual screen. All screens contain a wide 
range of turbulent structures, and these screens are carried 
by wind in different directions over the course of the ex- 
posure time. The distributions of wind speeds and atmo- 
sphere structure function parameters that PhoSim uses are 
based on observed data taken near the LSST site at Cerro 
Pachon, Chile. In P12, we demonstrated that with this at- 
mospheric model, the simulations from PhoSim show suf- 
ficiently realistic atmosphere-induced PSF spatial variation 
for the purpose of weak lensing studies. Although PhoSim 
was designed to simulate images from LSST in particular, 
the atmospheric model in PhoSim generates PSF patterns 
qualitatively generic to most large aperture telescopes. The 
result of this study can thus be easily extended to estimate 
the performance of our PSF interpolation method on other 
instruments. 

Finally, in this paper, we follow the convention in C12 
and simulate only the best 50% of the images in terms of 
image quality, where the median PSF size is ~ 0.7" full- 
width-half-maximum (FWHM). This is based on the fact 
that in previous cosmic shear measurements, most of the 
cosmic shear information appear to come from these "good" 
images ( jHoekstra et al.|2006[ ). That is, although one may be 
using all of the images, the bad images will be weighted low 
and thus the errors made in the PSF interpolation will also 
not be contributing as much to the final results. The range of 
atmospheric parameters used in the this work is consistent 
with C12, which also allows us to make statements of the 
spurious shear correlation function in Section[6.1| 



4.2 Testing programme 

We use a mock stellar catalogue to generate realistic images 
of star fields. The catalogue is based on the model of |Juric| 
et al. ( 20081), and contains a realistic population of stars in 



a typical LSST field with corresponding characteristics for 
each star. The average observed density of the population 
is ~ 1/arcmin 2 , which corresponds to that expected at a 
galactic latitude of \b\ ~ 60; the equatorial coordinates of the 
portion of the star catalogue used for this baseline simulated 
field was (1.5, +0.2) degrees. 

Two competing factors come into play in the PSF inter- 
polation problem: (1) the complexity of the PSF patterns, 
and (2) the number of stars available for constructing a PSF 
model (or effectively, the galactic latitude). The more com- 
plex the PSF pattern, or the fewer stars available to inter- 
polate, the more challenging it is to infer the PSF model 
from stars. The two effects are tested separately with the 
simulation and analysis pipeline described below. 

We address (1) by generating 100 realisations of the 
atmospheric PSF patterns, and "observing" the mock star 
fields at these 100 different "epochs." We generate these 
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Figure 3. Setting the entropic prior PDF on the map pixels. We show here in blue the entropic prior PDF (Equation [2f for a single 
pixel value h for two angular scales, 16 pixels (left) and 4 pixels (right). The grey curves are the pixel histograms of simulated ellipticity 
component maps at that scale, and the black curve is the average of those histograms. The entropic prior peaks at zero as required, 
and has slightly longer tails than a Gaussian. While the atmospheric PSF ellipticity map histograms have even longer tails and sharper 
peaks, especially on the smaller scales, by choosing the prior hyperparameter m judiciously, the width of the histogram can be matched 
reasonably well. Note that the two plots have different scales on the x-axis - the histogram is much tighter for the smaller scales. 



PSF realisations using the realistic distribution of atmo- 
spheric parameters described in P12. The median PSF size 
in our simulations is ~ 0.66" - for LSST, this corresponds 
to approximately half of the exposures with the best im- 
age quality. Each atmosphere realisation creates a unique 
PSF pattern; simulating 100 independently corresponds to 
a low-cadence survey campaign where observations are well- 
separated in time. 

To investigate (2) , we do not actually simulate star fields 
at different galactic latitudes. Instead, we create an over- 
sampled stellar catalogue based on the same stellar pop- 
ulation used for (1), so that our PhoSim input catalogue 
contains higher stellar density than an "average" field, while 
retaining the same signal to noise threshold in all "detected" 
star catalogs. We then down-sample at the detection cat- 
alogue level to achieve any desired stellar density used for 
interpolation - which can then be associated with a given 
galactic latitude. In this analysis, we consider stellar densi- 
ties between the range 4 and 0.25 arcmin -2 , which approxi- 
mately covers the range of galactic latitudes \b\ > 25. As will 
be explained later, the stellar density quoted here is after a 
signal-to-noise ratio (SNR) cut, which eliminates very dim 
stars and noise peaks (r >23.5). In reality, harder cuts may 
be used to guarantee purity of the star sample. 

For each of the above scenarios we generate one im- 
age containing only the expected stars at their given po- 
sitions, and a second image containing an ultra-dense grid 
of bright stars (~ 50/arcmin 2 ) to sample the "true" PSF 
pattern. Noise corresponding to a sky background level of 
22 mag/arcsec 2 was added to the star field images. All im- 
ages were generated in the r band for a single CCD sensor 
near the center of the LSST focal plane, which corresponds 
to a 13.6' x 13.6' field on the sky. 

The suite of simulated images were analysed using the 
same pipeline, in which stars were first detected using the 
Source Extractor ( Bertin & Arnouts 1996 1 package, and 



Kaiser]^] Shape estimation was performed using the imcat 
routine "getshapes". We use the output "e[0j" and "e[l]" as 
our representative measures of ellipticity (£i b fc and e^fe m 
Section^ and retained stars measured with imcat param- 
eter 'V larger than 25 (equivalent to a signal-to-noise cut 
~13). In the absence of uncertainty estimates on the ellip- 
ticity components, we propagate v as the stars' statistical 
weight. 

psfent, together with two other PSF modelling meth- 
ods, polynomial fitting and boxcar smoothing (details of our 
implementation of which are given in Appendix [A]| , were 
then applied to the detected stars. From each method, the 
output was a list of predicted values of the PSF elliptic- 
ity at the ultra-dense grid positions of the bright stars in 
the "true PSF" image. These grid positions stand for back- 
ground galaxy positions in a real lensing analysis. In this 
way we were able to compare the output catalogues directly 
with the true underlying ellipticity maps. 



4.3 Performance metrics 

To quantitatively evaluate the performance of the different 
PSF interpolation techniques, we employ two performance 
metrics in this paper. The first metric, <t[£psf], is defined to 
be the root-mean-square of the PSF ellipticity model error: 



^psf] = + (54) 



where 



Set 



— model £i,truc 



(5) 



(6) 



The second metric, 5" sys PS F, is defined as the average 
amplitude of the two-point correlation function of the PSF 
ellipticity model errors in the scales of interest: 



•'sya.PSF — 



\Z S + PSF (e)\de 



(7) 



then catalogued using the imcat software developed by Nick 



http : //www. if a.hawaii . edu/~kaiser/imcat/ 
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where £|_' PSF (#) is the correlation function of Se — Sei+i5e2- 
The absolute value in the integrand prevents the anti- 
correlation regime canceling out some of the correlation sig- 
nal. We use m in = 0.5' and # max = 13' in our single LSST 
CCD sensor (13.6' x 13.6') simulations. This range is cho- 
sen to sample the correlated PSF model errors on the full 
sensor while eliminating small scales far below the average 
stellar separation and large scales that come close to the 
boundaries. 

Note that ct[£psf] measures the level of the absolute 
errors in a certain PSF model, while <5"s ySi psp is a measure of 
the spatial correlation of the PSF model errors, in addition 
to their absolute levels. This motivates us to take the ratio of 
these two contributions and define an auxiliary metric F sys : 



°sys,PSF 

0"[£psf] 2 



(8) 



F sys is a relative measure of the how spatially-correlated the 
residual ellipticities are, independent of the absolute mag- 
nitude of the residual ellipticities. As seen in later sections, 
this figure facilitates our comparison of different PSF inter- 
polation methods. 



4.4 Scaling with number of exposures 

In the main analysis in this paper (Section [5j, we quantify 
the errors on the PSF ellipticity model for different interpo- 
lation techniques for a single LSST exposure. In reality, one 
can suppress the systematic errors that are independent be- 
tween frames when properly combining multiple exposures. 
It is the final combined systematic error of all the data that 
impacts the cosmological constraints from weak lensing. 

When the multiple exposures are far separated in time, 
the atmospheric conditions are different and one expects the 
PSF patterns to be independent. Similarly, the errors on the 
PSF model are also expected to be independent. This sug- 
gests that when averaging the ellipticity measurement of the 
same galaxy over N cxp exposures, the ellipticity errors from 
the atmosphere are expected to reduced as 1/ »/ N cxp (cor- 
responding to (t[epsf])j and the correlation of these errors 
should drop by l/JV ex p (corresponding to of. 



5"sys,PSF)i 



r[epgp] 



1 



(9) 
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<t[£psf] OC 



(10) 



This can be nicely demonstrated by performing the fol- 
lowing test: we take the 100 realisations of simulated atmo- 
spheric PSF and average the PSF model prediction at each 
position over the first jV oxp frames. o"[£psf] and ^ yB PSF is 
then calculated for the "average frame" and plotted against 
iV CX p in Figure [4] The two statistics scale with -/V eX p as ex- 
pected in Equations |9| and 1 1 0[ which confirms that the errors 
produced by all three interpolation methods used in this pa- 
per are stochastic over different exposures. 

We use the results in Figure [4] to support arguments 
later in Section [6. 1| where we estimate the number of expo- 
sures needed to reach a certain level of PSF model accuracy. 



Figure 6. Absolute two-point correlation functions of the PSF 
ellipticity model errors for the three different PSF patterns in 
Figure [5] In each panel, we show results for psfent (red), 5th- 
order polynomial fitting (green) and 5 X 5-pixel boxcar smoothing 
(blue). The hollowed labels indicate negative values. 



5 RESULTS 

5.1 Variation with PSF pattern 

The 100 atmospheric realisations in our simulation suite pro- 
vide a wide range of PSF patterns. For example, three very 
different PSF patterns on a single LSST CCD sensor are 
shown in the top row of Figure [5] where the colours repre- 
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Figure 4. (a) cfepgp] and (b) <T^ ys pgF calculated as a function of the number of exposures combined (averaged), -/V cxp . We show in 
both panels results for psfent (red), 5th-order polynomial fitting (green) and 5x5 boxcar smoothing (blue). The grey dash line indicates 
the l/y/ -/Vcxp and l/-/V eX p slope with arbitrary normalisation for (a) and (b) respectively. All three statistics in both plots roughly follow 
the iVexp scaling suggested by the grey lines. 



sent true PSF ei values. We treat the ei and £2 patterns as 
they were independent, as noted earlier, and only show ei 
here. 

We observe that the left column (a) contains a PSF 
pattern with some large-scale stripes that vary smoothly 
across the CCD sensor, with very fine ripples aligned at a 
direction different from the large stripes; the middle col- 
umn (b) contains medium-size blobs without any preferred 
direction; finally, the right column (c) shows nearly equal 
strength of stripes in two nearly orthogonal directions, cre- 
ating a grainy high-frequency pattern. When attempting to 
model these patterns at a typical galactic latitude of |6| ~ 60 
(stellar density ~ 1/arcmin 2 ), the stars available to us for 
reconstructing these PSF patterns are shown in the second 
row on the same ellipticity colour scale. The last three rows 
of Figure [5] show the PSF model generated from three dif- 
ferent PSF interpolation techniques in the following order: 
psfent, 5th-order polynomial fitting and 5 x 5-pixel boxcar 
smoothing. 

Visually, one can readily see the power of psfent in 
modelling the very different PSF patterns over the other 
two methods. In (a), all three methods failed to model the 
fine ripples, since they are much finer than the average stel- 
lar separations. They all do, however, manage to pick up the 
smooth "stripe" component. Both the polynomial and box- 
car model in this case show bad behaviours on the edges of 
the field, due to the small number of ill-measured stars domi- 
nating the model. In (b), where the underlying PSF patterns 
consist of mainly medium scale "patches", the Gaussian ICF 
used in psfent enables the model to capture these features 
nicely, which is not possible with a polynomial model. The 
boxcar model, on the other hand, is limited by the size of the 
filter, which in this case is slightly larger than the patches 
in the patterns. Finally in the right column (c), we show an 
example where the polynomial model becomes worse than 
even a simple boxcar smoothing. In this case the PSF pat- 



tern almost has no power on the large scales, causing the 
polynomial model to be entirely dominated by the noise. 

In all the cases shown here, we can see the flexible multi- 
scale algorithm allows psfent to model structures on a large 
number of scales, and is well regulated by the prior construc- 
tion thus less sensitive to noise. 

The absolute correlation functions |£1' PSF | for the three 
examples in Figure [5] are shown in Figure [6] In general, the 
model errors are more correlated on small scales due to the 
sparse sampling, and the limited resolution for all three mod- 
elling techniques. The shape of these correlation functions 
can be either smooth or oscillating. In particular, for poly- 
nomial models, the shape of £J PSF has some characteristic 
features: a transition from positive (correlation) to negative 
(anti-correlation) always appear at 3' - 4'. As discussed in 
H12, this is a result of both the modelling method and the 
true atmospheric PSF pattern. 

Figure [7] shows, for the 100 different realisations of the 
atmosphere, the median behaviour of these ellipticity error 
correlation functions with the error bars indicating the stan- 
dard deviation of the 100 curves divided by vT00. The me- 
dian o"[epsf] and cfsys.PSF values for these 100 atmospheric 
realisations are listed in Table [I] The final column Fays, as 
explained earlier, is a measure of the relative reduction in 
spatially-correlated residual ellipticity for the level of spatial 
correlation in the model errors independent of the absolute 
errors. Quantitatively, psfent provides ~ 17% improvement 
in the absolute ellipticity modelling error, or ct[£psf], over 
the 5th-order polynomial and ~ 22% improvement over the 
boxcar smoothing model. For the correlation of these er- 
rors, or 5"s ySi psF, psfent performs ~ 3.5 times better than 
the polynomial model and ~ 7 times better than boxcar 
smoothing. The corresponding F sys values suggest that the 
model errors from psfent is ~ 2.5 times less correlated than 
that from polynomial models and ~ 4.2 times less correlated 
than that from boxcar models. Notice that the main power 
in psfent lies not in the absolute reduction of the model 
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Figure 5. Illustration of the short exposure PSF interpolation problem, and the performance of different interpolation methods when 
we have very different PSF patterns. The three different realisations all have stellar density ~ 1/arcmin 2 . The maps in the first row 
show the "true" PSF ellipticity (ei) field that we would like to reconstruct from the stellar data in the second row, the observed stellar 
ellipticities. The last three rows show model PSF ellipticity maps constructed with psfent, a 5th-order polynomial fit and a 5x 5-pixel 
boxcar smoothing, respectively. 



© 2011 RAS, MNRAS 000,[T|fT5l 



Atmospheric PSF Interpolation for Weak Lensing 11 




O(arcmin) 



Figure 7. Median correlation function of the PSF ellipticity 
model errors of 100 different PSF patterns, at the "typical" stellar 
density of 1/ arcmin 2 . The results are shown for psfent (red) , 5th- 
order polynomial fitting (green) and 5x5 boxcar filtering (blue). 
The error bars indicate the rms spread in the 100 exposures di- 
vided by vlOO. Hollowed labels indicate negative values. 



sumption about the PSF pattern while ignoring the data. 
In contrast, the simple boxcar smoothing technique works 
in the opposite direction, where the model is purely driven 
by data with essentially no assumption on the expected PSF 
patterns. As a result, the models are just reflecting the avail- 
able data, where we get a model with no structure in the 
under-dense case (a), and a model with lots of small scale 
structure in the over-dense case (c). psfent is effectively a 
more sophisticated version of the boxcar smoothing with in- 
formative priors and multiple structure scales. The change 
from (a) to (c) for the psfent model is qualitatively similar 
to the boxcar smoothing, as it is primarily dictated by data. 
But when data is insufficient, as in (a), psfent does not 
generate entirely flat models like boxcar smoothing, rather, 
we can see traces of the psfent priors creating some struc- 
ture in the PSF model. When the stellar data is abundant 
(c), psfent lets data take over to drive the fit and only 
makes sure that the model stays in a physically reasonable 
range as specified by the priors. Figure [9] confirms the above 
observation more quantitatively. 



°1£PSf] °"sys,PSF F sys 

psfent 7.41xl0~ 3 2.41xl0~ 6 4.39xl0~ 2 

Polynomial 8.93 xlO" 3 8.67xl0~ 6 10.87xl0~ 2 

Boxcar 9.51 xlO -3 16.77xl0" 6 18.54xl0" 2 



Table 1. Median metric ct[£psf] an d cr 2 ys PSF for the three PSF 
interpolation techniques for 100 different PSF patterns sampled 
at the nominal stellar density of 1/arcmin 2 . The final column 
F syB is a measure of the level of the spatial correlations in the 
PSF model errors, independent of the absolute errors. 

errors, but in the fact that the flexible free-form model cre- 
ates makes errors less correlated in space, which is an impor- 
tant property for measurements like cosmic shear, where the 
main signal is embedded in the spatial correlation of galaxy 
shapes. 

5.2 Variation with stellar density 

Having examined the performance of the three PSF inter- 
polation methods on a "typical" field, we would now like to 
understand how the three PSF interpolation schemes are af- 
fected by the available density of stars (we explore the range 
from 0.25 to 4 stars per arcmin 2 for a complete sample of 
realistic stellar distributions). This test is especially impor- 
tant for images at high galactic latitude, where stars are 
very sparse. The ability to reconstruct the PSF variation 
at these fields may increase the effective area of a survey 
and therefore its statistical power. Figure [8] shows, for the 
PSF pattern in Figure [5](b), how the PSF model improves, 
for the three interpolation methods, as the stellar density 
increases. Figure [9] shows how the residual ellipticity corre- 
lation in each case changes accordingly. 

Figure [8] visually illustrates one example of how the 
three different interpolation methods respond to the in- 
creased available stellar data points. We observe that the 
polynomial models appear to be particularly ill behaved 
when the available stars are under-dense (a) and over-dense 
(c). This is an example of imposing an improper prior as- 



In Figure [To} we show for our 100 different atmosphere 
realisations the median <t[£psf] and 5f ys ,pgF statistics as a 
function of stellar density. In all stellar densities, psfent 
consistently performs ~ 20% better in ct[£psf] and 3-10 
times better in <5"s ySj psF compared to boxcar smoothing and 
polynomial fitting. 



The two statistics show similar trends in general. As we 
explained earlier, for psfent and boxcar smoothing, since 
the model is primarily driven by data, the model improves 
monotonically as more data becomes available. The 20% 
and nearly 4-times improvement of psfent in <t[£psf] and 
(5"sys,PSF respectively compared to boxcar smoothing mainly 
comes from psfent's ability to capture multi-scale struc- 
tures and regulate the model using priors so that noise does 
not get amplified. A fixed order polynomial function, on 
the other hand, when optimised for a certain stellar density 
(1/arcmin 2 ), over-fits (Figure [8] (c)) or under- fits (Figure [8] 
(a)) data when the stellar density varies, which results in 
a local minimum in the two green curves. The polynomial 
model, when optimised, is a reasonably good description of 
the smooth variation in the large-scale PSF patterns, but 
still fails to capture the small-scale structures, which ex- 
plains why psfent still performs better in that case. 

Although we emphasise the improvement of psfent 
over the other two methods, it is important to note at this 
point that the main improvement here is not from the spe- 
cific type of functional form (or un-parametrised model) one 
uses, but rather, the use of either realistic simulations or cal- 
ibration data to inform the prior PDFs for the flexible model 
parameters. We chose pixelated maps for use in psfent be- 
cause we learned from simulations that the PSF patterns are 
complicated and requiring a very flexible model. One can 
imagine, for example, an alternative method with the same 
spirit, where a basis of high-order polynomials are used to 
reconstruct the PSF variations with coefficients constrained 
by priors derived from simulations. 
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Figure 8. Illustration of the single short exposure PSF interpolation problem, and the performance of different interpolation methods as 
a function of stellar density. See the caption of Figure[5]for the description of the maps in the five rows; the different columns correspond 
to the same atmospheric PSF pattern (Figure [5] (b) ) , but sampled by stars with densities of 0.25 (a), 1 (b) and 4 (c) /arcmin 2 , as shown 
in the second row. 
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Figure 10. Median values for 100 PSF patterns for the two statistics: (a) o"[e P g p ] and (b) pSF over a range of stellar densities. The 
error bars show the interquartile range divided by \/100. We show in each panel results for psfent (red), 5th-order polynomial fitting 
(green) and 5x5 boxcar smoothing (blue). In both statistics, psfent performs consistently better than the other two techniques. 



6 DISCUSSION 



6.1 Implications of ct[£ P sf] and a 
systematics 



=.2 



PSP on shear 



To project the improvement on PSF interpolation from ps- 
fent onto the improvement in weak lensing shear measure- 
ments for a LSST-like survey, we return to Table [T] and dis- 
cuss the implications of the ct[£psf] and of ys PSF values un- 
der the nominal stellar density of 1/arcmin 2 . 



According to Amara & Refregier (2008 hereafter 
AR08) , for future stage IV ground-based weak lensing sur- 
vey^ jAlbrecht et al. 2006 1 not to be systematics-limited, 
one can set limits on the allowed systematic errors on the 
spurious shear power spectrum. [Paulin-Hcnriks son et aIT| 
( 2008 1 extended from AR08 and estimated that the allowed 



errors on determining the PSF ellipticity corresponding to 
those limits on spurious shear power spectrum is: 

-3 



o-[£psf] < 10" 



(11) 



Combining the first column in Table [T] and Equation |9j 
we can derive the number of exposures needed for each of 
the interpolation methods to achieve Equation [Tl] as listed 
in Table [2] Also listed is the corresponding operation time 
for LSST, where we have adapted the assumptions in C12 
and assumed two different scenarios - "optimistic" (A oxp = 
368) and "pessimistic" (N cxp — 184), where a total of iV cxp 
single 15-second exposures are combined in the final 10-year 
dataset for cosmic shear measurements. 

However, we note that Equation[TT]is a rather simplistic 
estimation that does not account properly for the correlation 
properties of these errors and the spurious shear arising from 
a realistic PSF correction pipeline. An alternative and more 
realistic approach to interpret the results from our study 
in terms of shear measurements is to turn to C12, where we 



9 In |Albrecht et al.| < |2006} , LSST is classified as a Stage IV 
"Large Survey Telescope (LST)" project. The Stage IV LST model 
assumes a ground-based survey with half-sky coverage, median 
readshift 1.0 — 1.2, and 30 - 40 well-measured weak lensing galax- 





iVexp 


operation time (years) 






optimistic 


pessimistic 


PSFENT 


55 


1.50 


2.99 


Polynomial 


80 


2.17 


4.35 


Boxcar 


90 


2.45 


4.89 



Table 2. Number of exposures required for the PSF ellipticity 
measurement accuracy to meet Equation |11| for the three PSF 
interpolation methods under nominal conditions. Also listed are 
the corresponding expected time span these exposures can be 
obtained by LSST in the optimistic and pessimistic scenarios. 





Nexp 


operation time (years) 






optimistic 


pessimistic 


PSFENT 


105 


2.86 


5.72 


Polynomial 


368 


10.0 


20.0 


Boxcar 


710 


19.34 


38.7 



les per arcmin 



Table 3. Number of exposures required for the PSF ellipticity 
measurement accuracy to meet the target value set by AR08 
and C12 for the three PSF interpolation methods under nomi- 
nal conditions. Also listed are the corresponding expected time 
span these exposures can be obtained by LSST in the optimistic 
and pessimistic scenarios. 



looked at the additive spurious shear correlation function by 
actually measuring these levels on high fidelity simulations. 
We concluded in C12 that the polynomial PSF model gen- 
erates spurious shear correlation approximately at the level 
required by AR08 in the optimistic case and 2 times too 
high in the pessimistic case. Since psfent provides a ~ 3.5 
times improvement in the PSF error correlation over poly- 
nomial models (Table [I] second column), we can expect it to 
also lower the spurious shear power spectrum by a factor of 
~ 3.5 if all the spurious shear is due to PSF interpolation er- 
rors. This brings the level of spurious shear power spectrum 
3.5 (optimistic) and 1.75 (pessimistic) times lower than the 
target level. Combined with Equation |10| we summarise in 
Table [3] the results in terms of the number of exposures and 
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survey time needed to achieve the target value set in AR08. 



7 CONCLUSIONS 



On the other hand, if one considers that there are shear 
measurement algorithm errors that have not been accounted 
for in C12, then the situation could be worse. Assume, for 
example, these unaccounted algorithm errors take up to 50% 
the allowed systematic errors set by AR08, the improvement 
in psfent then results in a spurious shear power spectrum 
about 1.3 times lower (optimistic) and 1.1 times higher (pes- 
simistic) than the target level. In that case, even psfent will 
be a marginal failure in the pessimistic scenario. 



6.2 Correlating galaxies across exposures 

Correlating galaxy ellipticities in different exposures has 
been one of the proposed solutions to the problem of stochas- 



tic atmospheric PSF correlations (Jain et al. 2006). This, 



however, comes at the price of decreasing the statistical 
power of the survey. Exploring the tradeoff between statisti- 
cal and systematic PSF interpolation errors using this tech- 
nique is not the main focus of this paper. Here, we have effec- 
tively assumed a LENsFlT-style ( |Miller et al.|2007| ) analysis, 
where the PSF ellipticity map is estimated for each expo- 
sure, the galaxy images deconvolved, and the ellipticity es- 
timates simply averaged. This assumption is justified by the 
behaviour of the residual PSF ellipticity correlation function 
with increasing N cxp as we have shown. We presume that 
when correlating between different exposures, it will still be 
beneficial to start with a more accurately-interpolated PSF 
model. 



6.3 Computational cost 

On a standard 64-bit 2x2.2 GHz processor, psfent takes 
on average ~13.0 seconds per LSST CCD sensor per pair 
of shape parameters (£1,62). Adding error estimation (sam- 
pling from the posterior cloud) adds an extra ~9 sec of 
run time. This is about an order of magnitude slower than 
the two other techniques we investigated. For more com- 
plicated shape parametrisation the runtime would increase 
linearly with the number of parameters, and also with the 
number of images analysed. For example, if psfent were to 
be used for in the LSST data processing pipeline, the PSF 
model in each exposure on the ~200 CCD sensors would 
need to be calculated in ~ 19 seconds (15-second exposure 
+ 2-second readout + 2-second shutter open/close). This 
demands ~200 computers running for the PSF reconstruc- 
tion alone. Though large, these numbers are not outrageous 
considering the expected decrease in unit cost for computers 
over the next decade and the possibility of further acceler- 
ating the code via hardware parallelisation such as CPUs. 

In the meantime, we recommend that the current code 
can be used for smaller scale datasets when extremely ac- 
curately interpolated PSF maps are vital for the specific 
science goal. As well as weak lensing, these might include 
reconstructing the detailed structures of complex objects, 
and precision photometry of faint objects. 



In this paper we have introduced a new PSF interpolation 
method, psfent, based on a multi-scale maximum-entropy 
image reconstruction code. The problem we set out to solve 
is reconstructing the multi-scale spatial variations of PSF 
shapes, due to atmospheric turbulence, in short exposure 
images, from sparsely distributed, noisily measured stars - 
a potential problem for future weak lensing surveys such as 
LSST. 

Our analysis of simulated data allows us to draw the 
following conclusions: 

• Compared with two other PSF interpolation methods, 
one of which is commonly used in current data analyses, 
psfent provides more accurately interpolated PSF shapes: 
the absolute residual ellipticities improve ~ 20% while the 
correlated residual ellipticities are a factor 3 -10 smaller than 
the other two methods, over a wide range of different PSF 
patterns and stellar densities. 

• When combining multi-epoch datasets, the interpola- 
tion errors due to the atmosphere are stochastic, decreasing 
with the number of exposures taken as 1/ ^/JV CX p. The cor- 
relation function amplitude decreases as 1/N exp . 

• The improvement in PSF modelling from psfent sup- 
presses the spurious shear in cosmic shear measurements. 
Combining previous studies from realistic simulations and 
the results in this work, we estimate that systematic errors 
due to PSF interpolation for LSST will be 3.5 (optimistic) 
and 1.75 (pessimistic) times lower than the statistical errors. 

• Taking into account other algorithm errors, however, 
the improvement from psfent may not be sufficient to bring 
the systematic errors down to the target level in the most 
pessimistic scenario. 

While it may still become a practical solution even for 
large datasets such as LSST, the relatively high computa- 
tional cost of powerful algorithms like psfent should moti- 
vate the development of faster inference methods. The paper 



by Berge et al. (2012 I appeared as we were completing this 



study: their results are complementary to those presented 
here, in that they investigate the PSF interpolation problem 
on larger scales, using re-sampled real Subaru data that is 
less dominated by the atmosphere. It would be very interest- 
ing to apply their Kriging technique to our simulated data, 
and compare performance in terms of accuracy and speed. 
Like |F3erge et al.] p012| , we have demonstrated the benefits 
of using very flexible models, but also the use of methods 
both inspired and constrained by realistic simulations. In- 
jecting statistical information about the atmospheric PSF 
anisotropy into interpolation methods appears to be a fruit- 
ful approach. 

All the simulated data used in this paper is freely avail- 
able at the following website: http : //www. slac . Stanford . 
edu/ ~ chihway /psf ent / . 
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APPENDIX A: OTHER MODELLING 
METHODS 

Al Polynomial fitting 

As discussed in Section [2] in most current weak lens- 
ing pipelines, the PSF spatial variation is assumed to be 
smoothly varying on scaled comparable to the field and can 
thus be modelled with some low order polynomial functions. 
We have shown that PSF patterns from short exposures, 
however, display higher frequency spatial variations that re- 
quire higher order fits. As a result, we have examined polyno- 
mials of order 2 to 5 as PSF models and found that 2nd-order 
models are insufficient in representing the PSF variations; 
3rd-order models are usually sufficient, but occasionally an 
even higher order model is needed. We choose to use 5th- 
order models in our main analyses, but have confirmed that 
in most cases, the results are identical to 3rd-order fits. We 
minimise the effective x 2; 

X^E^r^P;^)-^) 2 , i = l,2, (Al) 

where e™ odel ' : ' {p;x,y) is a two dimensional, jth order poly- 
nomial function of (x,y), p are the fitting parameters and 
we use the signal-to-noise ratio of each star as w E . 

A2 Boxcar smoothing 

Alternatively, we examine the approach of modelling the 
PSF by directly smoothing the stellar ellipticities with a 
boxcar filter. In this approach, stars need to first be binned 
into large pixels, where we assign a weighted average ellip- 
ticity to each pixel (j\,]2), defined: 

sp= = Y, WeSi ' * = 1 > 2 - 

pixel (ii,j 2 ) 
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The size of the pixels are determined by the number of stars 
- we pixilate the image so that in each pixel contains ap- 
proximately 1 star. A boxcar filter of size m x m pixels is 
then applied to the large pixel grid so that the ellipticity 
of pixel is replaced by the average ellipticity of the 

neighbouring (m 2 — 1) pixels: 



_model,j 1 j 2 



tr« E E 4 



i m-1 . 
m — — - — , i = l,2 



The original coordinates are then recovered so that the mod- 
elled PSF ellipticities of all points falling in pixel (ji, J2) will 
be £ model 'iii2 



APPENDIX B: PSFENT WITH OTHER SHAPE 
PARAMETERS 

It had been shown in H12 that due to atmospheric distor- 
tion, the size of the PSF varies on spatial frequencies similar 
to that of £1 and £2. This suggests that with the same pro- 
cedure described in Section [4] we could use psfent to inter- 
polate the PSF size, and by extension other shape parame- 
ters. We demonstrate below our results of interpolating the 
shape parameter "FWHM WORLD" output from Source 
Extractor (which is an estimate of the PSF FWHM size in 
arcseconds) with the three interpolation methods used in 
Section [5] 

First we derive the appropriate priors (to* values) 
for FWHM _ WORLD from the "true" PSF maps. We 
have m i= [0.0001, 0.0002, 0.0026, 0.0055, 0.0079, 0.0100, 
0.0206] for FWHM WORLD. The set of priors for the 
FWHM WORLD variation appear to follow the same trend 
as that of £1 and £2, while putting slightly more power 
on the larger scales. This can also be observed qualita- 
tively in Figure |B1| where we show the equivalent of Fig- 
ure [5] (b), with the interpolant being the PSF size, R (or 
FWHM _ WORLD), instead of £1. We then calculate for 
the 100 atmosphere realisations and for the three different 
PSF interpolation methods, the statistics in analogy with 
<t[epsp]: 



a[8R 2 \ _ Vi(SR 2 ) 2 ) 



(R 2 ) 



(R 2 ) 



where 



8 R 2 — \R 2 no d e i 



— R tr i 



(Bl) 



(B2) 



We find that the median A values for the 100 simulations 
are 1.9xl0 -2 for psfent, 2.9x10 -2 for 5th-order polyno- 
mial models and 4.7xl0 -2 for 5x5 boxcar smoothing. That 
is, psfent provides an improvement of ~ 35% over polyno- 
mials and ~ 60% over boxcar smoothing in interpolating the 
PSF size. Note that the A value for polynomial fitting ap- 
pears much larger than that derived in H12. This is because 
the atmospheric variations in our 100 images is much larger 
than the data used in H12, where they have used 33 images, 
all taken within 2 nights. We conclude here that psfent 
is capable of interpolating other general shape parameters, 
once we properly informed the priors from simulations. 



APPENDIX C: CORRELATION BETWEEN 
ELLIPTICITY COMPONENTS 



We have argued in Section |37lj that the two ellipticity compo- 
nents should be independent from each other in each expo- 
sure, while both varying with similar amplitudes and spatial 
structures between different exposures. In this appendix we 
demonstrate this by calculating the following three spatial 
correlation functions for the 100 "true" PSF images: 



£u(0) = <ei(0o)ei(0o + 0)); 

e22(0)=<£2(0 O )£2(0„+0)>; 

£i 2 (0) = (ei(0o)ea(0o + 0)>. 



(CI) 



(C2) 
(C3) 



In Figure [CT| the three correlation functions are shown, 
each representing the median of the 100 realisations in the 
sample, with the error bars showing the standard deviation 
of the 100 realisations divided by Vl00. The £n and £22 
curves are clearly positive and show similar structures. The 
difference in the two curves may be due to sample variance, 
and the fact that there are no rotation/dithering between 
these exposures. On the other hand, £12 is consistent with 
zero. This supporting our argument that the "phase" of the 
two ellipticity components are independent of each other, 
while the spatial power spectrum is similar. 
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Figure 9. Absolute two point correlation functions of the PSF 
ellipticity model errors for different stellar densities: (a) 0.25, (b) 
1 and (c) 4 /arcmin 2 . The three cases have the same underlying 
PSF pattern but are sampled at different rates. In each panel, we 
show results for psfent (red), 5th-order polynomial fitting green) 
and 5x5 boxcar smoothing (blue). The hollowed labels indicate 
negative values. 



r 
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Figure Bl. Illustration for interpolation of PSF size (FWHM 
size in arcseconds). The top map shows the "true" PSF size field 
that we would like to reconstruct from the stellar data in the 
second map, the observed stellar size. The next three maps show 
model maps constructed with psfent, a 5th-order polynomial fit 
and a 5x 5-pixel boxcar smoothing, respectively. The atmospheric 
PSF pattern and stellar density is the same as that in Figure [5] 
(b). 
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Figure CI. Three diagnostic correlation functions fn (red), £22 
(green) and £12 (blue), plotted with the full ellipticity correlation 
function §+ (grey). Each curve represents the median of the 100 
different atmosphere realisations. The error bars indicate the rms 
spread in the 100 exposures divided by V100. 
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